/*
 * Tracker.h
 *
 *  Created on: Nov 28, 2011
 *      Author: schurade
 */
#ifndef TRACKER_H_
#define TRACKER_H_

#include <string>
#include <vector>

#include "nifti/nifti1_io.h"

using namespace std;

class Tracker {
public:
	Tracker( size_t minLength, float minFA, float stepSize, float smoothness );
	virtual ~Tracker();

	bool load( vector<string> filenames );
	bool loadMask( string filename );

	void fitTensors();

	void saveTensors();
	void saveFA();
	void saveEvec();

	void writeNifti( string fn, vector<float>data, int dim );

	void saveFibs();

	void computeTrace();

	void computeFA();

	void computeEigen();

	void startTracking();

private:
	vector<float> track( size_t x, size_t y, size_t z, bool negDir );

	size_t getID( float x, float y, float z );

	vector<float> getInterpolatedTensor( size_t id, float x, float y, float z );

	void createLogTensors();

	vector<float>logT( size_t id );

	vector<float>expT( vector<float> t );

	// input file related
	vector<string>filenames;
	vector<float> bvecs;
	vector<int> bvals;
	nifti_image* header;
	nifti_image* maskHeader;
	double maxEval;
	float evalScale;

	// loaded values from nifti header
	size_t numGradients;
	size_t vDim;
	size_t numB0;
	size_t columns;
	size_t rows;
	size_t frames;
	float xVoxelSize;
	float yVoxelSize;
	float zVoxelSize;

	size_t blockSize;

	size_t m_minLength;
	float m_minFA;
	float m_stepSize;
	int maxStepsInVoxel;
	float m_smoothness;

	// hardii file split up
	vector<float> b0Images;
	vector<int16_t> dti;

	// calculated tensors, fa and eigen vectors
	vector<float>tensors;
	vector<float>trace;
	vector<float>fa;
	vector<float>eval1;
	vector<float>eval2;
	vector<float>eval3;
	vector<float>evec1;
	vector<float>evec2;
	vector<float>evec3;
	vector<unsigned char>m_mask;
	vector<float>logTensors;

	vector< vector< float > >fibs;

	double pow2( double v ) {return v*v;};
	double pow3( double v ) {return v*v*v;};
};

template< class T > T switchByteOrder( const T value )
{
    size_t numBytes = sizeof( T );
    T result = value;
    if( numBytes == 1 )
    {
        return result;
    }
    char *s  = reinterpret_cast< char* >( &result );
    for( size_t i = 0; i < numBytes / 2; ++i )
    {
        std::swap( s[i], s[ ( numBytes - 1 ) - i ] );
    }
    return result;
}

/**
 * Transform a whole array of elements (of type T and size of sizeof(T))
 * into opposite byte order.
 *
 * \param array Array containing the data
 * \param arraySize The number of elements which is not the number of
 * bytes but e.g. the number of floats
 */
template< class T > void switchByteOrderOfArray( T *array, const size_t arraySize )
{
    for( size_t i = 0; i < arraySize; ++i )
    {
        array[i] = switchByteOrder< T >( array[i] );
    }
}

#endif /* TRACKER_H_ */
